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Abstract 

Within the nonrelativistic QCD (NRQCD) factorization framework, we investigate the inclusive 
production of the h c meson associated with either light hadrons or charmed hadrons at B factory 
energy yfs = 10.58 GeV. Both the leading color-singlet and color-octet channels are included. For 
the h c production associated with light hadrons, the total production rate is dominated by the 
color-octet channel, thus the future measurement of this process may impose useful constraint 
on the value of the color-octet matrix element (Og c ( 1 S'o)); for the h c production associated with 
charmed hadrons, the total production rate is about one order of magnitude smaller, and dominated 
by the color-singlet channel. 
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I. INTRODUCTION 



The lowest-lying cc( 1 Pi) state, the h c (lP) meson, is the last one found among all the 
charmonium members below open charm threshold. This elusive particle was not firmly 
established until 2005, through the isospin-violating decay ip' — > 7r°(— > 77)/i c ( — > 7^7) by 
CLEO Collaboration [1], as well as through the process pp — >■ h c — >■ 77^,7 by the E835 
experiment [2]. Quite recently, the analogous 1 Pi members in the bottomonium family, the 
h b (lP, 2P) mesons, have also been observed by the Belle Collaboration through the process 
e+e- ->■ T(5S) ->■ h b (nP) + tt+tt- [3]. 

The quite accurate measurements of the mass of the h c (also h b (lP, 2P)) implies a rather 
small P-wave hyperfine mass splitting. This is theoretically intriguing since it might be 
able to impose some severe constraint on the quark spin-spin interaction as well as possible 
charm meson loop effect. 

Aside from its mass [4, 5], our knowledge about the h c state is still quite limited. The 
h c appears to be a narrow resonance with the total width of 0.73 ± 0.45 ± 0.28 MeV [5]. 
So far, only two decay channels of the h c have been measured, one is the hadronic decay 
h c — > 2(7T + 7r~)7r° [6], and the other is the much more abundant El transition h c — > r\ c ^. 
Recently BES III experiment has measured the absolute branching fraction of the latter 
process and gives B(h c — > 7^) = (54.3 ± 6.7 ± 5.2)% [5]. 

In contrast to the decay, our understanding of the h c production is even poorer. The 
only measurement of the h c production is from a recent CLEO experiment, by observing 
the process e + e~ — >■ h c n + n^ at y/s = 4.170 GeV [7]. On the theoretical side, rather few 
works on h c production are scattered in literature. Among these studies, are h c production 
in B meson inclusive decay [8, 9], h c photoproduction [10], h c hadroproduction [11, 12]. It 
is interesting to compare this situation with highly intensive studies on J/ip production in 
various collision experiments [13]. 

Our goal in this work is to carry out a detailed study on inclusive h c production in 
e + e~ annihilation, which is specifically relevant to the B factory experiments. Specifically, 
we investigate the inclusive production of the h c meson associated both with light hadrons 
and with charmed hadrons. There are both theoretical and experimental merits to study 
these processes. On the theoretical side, since these processes are much simpler than the 
h c production in hadronic collision, one expects to obtain more precise prediction with less 
contamination from the nonperturbative side of QCD. On the experimental side, since e + e~ 
collision experiment possesses much clean background than the hadronic collider, and the 
B factories have already accumulated a quite large data sample near the T(4S) resonance, 
it could well be the most likely place to unambiguously observe the elusive h c signal. 

Our study is based on the nonrelativistic QCD (NRQCD) factorization approach, a 
widely-accepted framework to deal with inclusive quarkonium production [14] , which heavily 
exploits the nonrelativistic nature of heavy quarkonium. A highlight of this approach is the 
so-called color-octet mechanism, which serves as an indispensable ingredient in order to give 
a meaningful prediction for the P wave quarkonium production [8]. 

As we will see, for the process e + e~ — > h c + light hadrons, the color-octet mechanism 
indeed plays a pivotal role in rendering infrared-finite prediction for the inclusive h c produc- 
tion rate. Furthermore, this process is found to be dominated by the color-octet channel. 
Our study reveals that this process may have a sizable total cross section that is comparable 
in magnitude with that of e + e~ — > J/ip+ light hadrons, which has been measured some 
time ago by the B factory experiments [15]. Future measurement of this process at the B 
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factories may provide an explicit test on the color-octet mechanism, as well as put some 
useful constraint on the value of the color-octet matrix element for h c 1 . 

The rest of the paper is organized as follows. In Sec. II, we present the NRQCD factoriza- 
tion formulas for h c inclusive production in e + e~ annihilation, accurate at lowest order in v 
(the characteristic velocity of charm quark inside the h c meson). In Sec. Ill, by employing the 
perturbative matching ansatz, we determine the infrared-finite color-singlet and color-octet 
short-distance coefficients associated with the process e + e~ — > h c + light hadrons. In Sec. IV, 
we determine the color-singlet and color-octet short-distance coefficients associated with the 
process e + e~ — > h c + charmed hadrons. Based on our calculation in previous sections, we 
devote Sec. V to exploring the observation prospects of h c production at B factories. Finally 
we summarize in Sec. VI. In Appendix A, we explain some technical details in isolating the 
infrared divergence for e + e~ — > h c + gg in dimensional regularization. In Appendix B, we 
present the color-singlet short- distance coefficient for the process e + e~ — > r\ c + cc. 



II. NRQCD FACTORIZATION FORMULA FOR h c PRODUCTION 

NRQCD factorization formalism is a systematic tool for analyzing the inclusive produc- 
tion of heavy quarkonium [14] . The production rate can be expressed as a sum of products 
of short-distance coefficients and nonperturbative, albeit universal vacuum NRQCD matrix 
elements, whose importance is organized by the typical quark velocity, v. In this work, we 
will consider e + e~ — > h c + X in this factorization framework. At the lowest order in v, the 
velocity counting rule implies that the cross section of h c has the following form: 

da[e + e - ^ hc + X } = —L^CP,)} + _!<Cfcf So)), (1) 

m c m c 

where the color-singlet operator C^ C ( 1 P 1 ) and the color-octet operator Og c ( 1 5'o) have been 
introduced in [14]. It is interesting to contrast this P-wave quarkonium production process 
with the «S-wave onium production, where the color-octet operator matrix element first 
comes into play only at relative 0{v A ). 

dFi, dF$ are infrared-finite short-distance coefficients associated with the vacuum matrix 
elements of color-singlet and -octet NRQCD production operators. Since they are insensitive 
to the long-distance strong interaction dynamics, a standard way of determining them is 
through the perturbative matching procedure: replacing the h c state appearing in (1) with 
the free on-shell cc(n) states (n = 1 S'q 8 ' ) or 1 P^), and computing both sides of (1) using 
perturbative QCD and perturbative NRQCD, respectively, then enforcing that they generate 
identical results. Finally, one then solves two linear equations to identify the two unknown 
coefficients, order by order in q s 2 . 

The perturbative matching for e + e~ — > h c + charmed hadrons is straightforward. In 
contrast, the matching procedure for e + e~ — > h c + light hadrons is more subtle and involved, 



1 The majority of the results in this paper has already been presented in Ref. [16]. 

2 In this work, we are only looking for the h c energy distribution and total cross section, rather than its 
angular distribution. To this purpose, we may adopt a standard shortcut to simplify the intermediate 
calculations [17]. First compute the virtual photon decay into h Cl then use the following formula to convert 
the decay rate into h c cross section: da[e+e~ -t h c (P) + X]/dP° = dT[<y* -> h c (P) + X]/dP°, where 
P M represents the 4-momentum of the h c in the e + e~ center-of-mass frame. 
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since the QCD side calculation for the color-singlet channel da[e + e~ — > cc( 1 P 1 ( ') + gg] 
contains infrared divergence, which must be absorbed into the color-octet matrix element 
to render an infrared (IR) finite color-singlet short-distance coefficient. 

III. h c PRODUCTION ASSOCIATED WITH LIGHT HADRONS 

In this section, we will apply the perturbative matching procedure to deduce the short- 
distance coefficients for inclusive h c production associated with light hadrons at B factories, 
i.e., e + e~ — > h c + light hadrons. For clarity, we will always attach a superscript "LH" to 
both Fi and F 8 in this section, to differentiate them from the analogous coefficients for the 
h c production associated with charmed hadrons, which will be reported in the next section. 

A. Determining dFg U 

We begin with calculating the differential short-distance coefficient <iF g LH affiliated with 
the color-octet operator. At the lowest order in a s , only two Feynman diagrams need to be 
considered for the process e + e _ — > cc^S^) + g. 

As mentioned before, the differential cross section da[e + e~ — > cc( 1 P 1 ) + gg]/dz would 
develop an IR divergence as one of the gluons gets soft. In order to deduce the IR finite 
color-singlet coefficient gLF\ lh , the perturbative factorization formula (1) implies that we 
should also consider the color-octet coefficient dF$ multiplied by the renormalized 0(a s ) 
color-octet matrix element, which is generally IR divergent. Throughout this work we find 
it most convenient to employ the dimensional regularization (DR) to regularize both UV 
and IR divergences. Therefore, it is necessary to compute e + e~ — > cc( 1 S^) + g in D = 4 — 2e 
spacetime dimensions. 

It is convenient to use the covariant spin projection method [18] to compute the amplitude 
of 7 * -> ccW) + 

g in D spacetime dimensions. A subtlety is that the appearance of 75 
from the spin-singlet projector, which requires some care to handle it in D dimensions. 
For consistency, we adopt the 't Hooft-Veltman (HV) prescription [19], and utilize West's 
formula to calculate the trace involving one 75 and a string of Dirac matrices [20]. The 
Levi-Civita tensor is assumed as a 4-dimensional object 3 . 

We will always work in the e + e~ center-of-mass frame throughout this paper, where \fs 
denotes the e + e~ center-of-mass energy. For notational simplicity, we define the energy 
fraction of h c , z = 2P°/ v /s, as well as the ratio of charm quark mass over center-of-mass 
energy, r = 4ml/ s. The differential two-body phase space in D dimensions can be expressed 
as 



3 As a crosscheck, we have also tried to treat the Levi-Civita tensor as a D-dimensional object when 
computing the squared amplitude. After matching is done, one is justified to return to 4 dimensions. It 
turns out that this alternative prescription leads to the identical short-distance coefficients dF% and dF\, 
as given in (4) and (17) in the text. 



d$ 2 = ^s~ e (l - rf- 2e 5{\ + r - z)dz, 



(2) 



where c t = (47r) 



e r(i-e) 

r(2-2e) • 
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Comparing both sides of (1) at 0(a s ), it is easy to find the differential coefficient dF 8 LH 
in D dimensions: 

rfF 8 LH / y 2 \ e 'Z2ii 2 e 2 c a 2 a s m u _ 2eA/1 , , 

~ c A — \ ^2 i 1 -^ S(l + r-z), (3) 



dz V s / 3s 2 

where \i is the compensating mass scale in DR, e c = | is the electric charge of charm quark. 
In the limit e — > 0, the differential color-octet coefficient reduces to 

(l--)*(l + r-.). (4) 

B. Determining cLF\ LH 

We proceed to determine the color-singlet coefficient for h c production associated with 
light hadrons in e + e~ annihilation, dF^. This can be expedited by replacing the h c state 
with a free cc( 1 P 1 <1 ' ) ) pair, and matching both sides of (I) that are computed in perturbative 
QCD and perturbative NRQCD, respectively 

At the lowest order in a s , the color-singlet h c production associated with light hadrons 
can proceed through the parton-level process e + e~ — > cc(}P^)gg, which can be picturized 
by six Feynman diagrams. Let P, ki, k 2 signify the momenta of the cc{}P^) pair, gluon 1, 
and gluon 2, respectively. We will always assume P 2 ~ Am 2 . It is convenient to introduce 
three fractional energy variables z, x\ and 

2P° _ 2A;? _ 2A;° 

z — — X\ — —=, x 2 — — p, (5) 
V s V s V s 

which are subject to the constraint x\ + x 2 + z = 2, as required by the energy conservation. 

Since we employ the DR to regularize the potential IR divergences, it is useful to write 
down the 3-body phase space integral in 4 — 2e spacetime dimensions: 

where 9 signifies the angle between P and k x in the 7* rest frame. For given z and x 1: it 
can be uniquely determined: 

2(1 +r -: ] - xi(2 - z) ( 7) 

x\\/ z 2 — 4r 

Since we have taken the shortcut by first calculating the virtual photon decay into 3- 
body final state, only two independent energy fraction variables need be retained in the 
integration measure in (6). The integration limits of z have been explicitly labeled in (6), 
while the integration boundaries of #i, xf, are parameterized as a(z) ± b(z), where 

a(z) = \{2-z), (8a) 
b(z) = \^z 2 - 47. (8b) 
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We again use the spin projection technique [18] to compute the amplitude of 7* — > 
cc( 1 P 1 ^ 1 ' ) ) + gg in D dimensions and use the HV prescription to handle the trace involving 
75. After squaring the amplitude and summing over polarizations and colors, we separate 
the squared amplitude into two pieces: 



A[y* ->■ ccCP^) + gg}~ = /div(si, z) + l &Q { Xl , z), 



Pol, Col 



(9) 



where 



+ 



(1 + r — 2 — X1) 2 (1 + r — z — x 2 ) 



(10) 



represents the term that would bring forth an IR singularity when integrating over two 
distinct phase space corners: z — > 1 + r, x 1 — >■ and 2 — >■ 1 + r, x 2 — > 0, where one of 
the gluons become soft. The symbol ifi n denotes the remainder of the squared amplitude 
which contains no terms as dangerous as those in (10), therefore renders a finite result 
upon integrating over the entire 3-body phase space. Bose symmetry guarantees that ig n is 
manifestly symmetric under interchange between x\ and xi — 1 — z — x\. Since its explicit 
expression is somewhat lengthy, so will not be reproduced here. 

Accordingly, the energy distribution of the parton cross section can also be decomposed 
into two parts: 



da[e + e- ->■ cc^P^, P) + gg] 
dz 



da div da &n 



dz 



dz 



(11) 



which are obtained by integrating (9) over the entire momentum range of gluon 1: 



da, 



rl+r 

J2^F dz 
-1+r 



div 



3s 



2 , d® 3 I diY (x 1: z), 



J2^7- dz 



^ J d^ 3 h n (x u z). 



(12a) 
(12b) 



In deriving these, we have used the conversion formula explained in footnote 2, as well as 
included a factor ^ to account for the indistinguishability of two gluons in the final state. 

It is straightforward to complete the integration over x\ in the right side of (12b). Since 
everything is finite, this integral can be directly calculated in 4 dimensions. In contrast, 
integrating J div over x\ in DR requires some special care due to emergence of the IR diver- 
gence that occurs at z — 1 + r. We devote Appendix A to expounding the intermediate 
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technical steps, and here just simply jump to the desired results: 
d&dn 1287T e 2 a 2 Cpa 2 c e (47r)' 



'div 

dz 



3m 2 s 2 



r(i-e) 



em 4m^ Vr 



1 - r)5(l + r — z) x 
1 



1 + r — z 



4r 



£ — 2r 



+ Vz 2 - 4r 



2z 



2 Z 



4r 



2; — 2r l + r — z 

d(7fi n 2567re 2 a: 2 CVa: 2 



(13a) 



1 



2r)Vz 2 — 4r x 



3m 2 s 2 (2 - z)\z - 2rf 

16r(-3 + 2r - 6r 2 - 6r 3 + 3r 4 + 6r 5 ) + 16r(7 + 2r + 24r 2 + r 3 - 20r 4 - 2r 5 )z 
+4(2 - 35r - 72r 2 - 74r 3 + 94r 4 + 25r 5 )z 2 - 8(2 - 21r - 34r 2 + 18r 3 + 15r 4 )^ 3 
+2(3 - 59r - 6r 2 + 32r 3 )/ + (3 + 23r - 14r 2 )^ 5 - z 6 
2r + Vz 2 - 4r 



+ ln- 



- 32r 2 (3 - r + 4r 2 + 3r 3 + 2r 5 + 5r b ) 



z — 2r — V-2 2 — 4r 
+32r 2 (8 + r + 14r 2 + 4r 3 + 6r 4 + 21r 5 + 2r 6 )z 
+8r(4 - 37r - 31r 2 - 54r 3 - 38r 4 - 149r 5 - 31r 6 )^ 2 
-8r(12 - 25r - 26r 2 - 42r 3 - 148r 4 - 51r 5 )^ 3 

+2r(61 - 17r - 55r 2 - 363r 3 - 186r 4 )^ 4 - 12r(8 + 2r - 21r 2 - 17r 3 )z 5 



+ (1 + 45r - 37r 2 - 65r 3 )z 6 - (1 + 5r - 10r 2 )z 7 



(13b) 



The appearance of the (^-function and "+" -function exhibits some peculiarity of the energy 
distribution da^jdz near the maximal h c energy. As usual, these functions should be 
interpreted as the distributions in the mathematical sense. The [/(z)] + function in (13a) 
is defined such that when convoluting it with an arbitrary function g(z) that is regular at 
z — l + r, one gets 

/•l+r rl+r 

/ dz[f(z)] + g(z) = dzf(z)(g(z)-g(l + r)). (14) 

J2jr J2y/r 

From (13a), we find that the IR singularity is exactly located at the maximal value of z. 
One certainly expects that this IR singularity will be swept out once including the color-octet 
contribution. Encouragingly, equation (3) implies that the differential color-octet coefficient 
dFg U is also proportional to a peaked distribution 5(1 + r — z). 

Since dF$ n /dz is of order a s only, in order to match the 0(a 2 ) accuracy of the color- 
singlet parton cross section, we need incorporate the 0(a s ) correction to the perturbative 
color-octet NRQCD matrix element. This correction has already been inferred previously 
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in the related work on quarkonium production, so we just present the result 4 : 



(OrCS ))us = (OfCS )) 



AO) 



2C F a s 
3N c Trm 2 c 



— + In An 

em 



(o) 



(15) 



where Cp = N ^ N 1 , and N c = 3 is the number of the colors in QCD. Under renormalization, 



2JV, 

the color-octet operator O^^Sq) mixes with the color-singlet operator (9J C ( 1 P 1 ) (it also 
mixes with the color-octet operator Cg c ( 1 Pi), but whose effect would arise at higher order in 
v). It is important to note that, the MS-renormalized perturbative matrix element (Og c ( 1 S'o)) 
develops a logarithmic IR divergence. 

In passing, it may be worth stressing that, unlike the NRQCD decay operators, the 
analogous production operators are no longer the local 4-fermion operators. In fact, in a 
series of work [21-23], Nayak, Qiu and Sterman have recently advocated that one must insert 
the proper gauge links to the original definitions of the color-octet production operators [14] 
to warrant the gauge invariance, and the nontrivial effect due to this gauge completion first 
shows up at next-to-next-to- leading order in a s (it has also been explicitly examined that, at 
the NLO in a s , forgoing the gauge completion for the color-octet operators does not bring in 
any inconsistency [24]). Thus to our purpose, we are content with staying with the original 
definition given in [14]. In this respect, equation (15) looks very similar to the analogous 
formula for the renormalized color-octet decay operator in NRQCD [18, 25]. 

In light of the matching condition in (1) for the cc(}P^) channel, one can write down 
the following equation for dF^ to the order a 2 : 



da[e + e- ->■ ccQP^) + gg\ _ 6JV c dF 1 LH (^) 



dz 

1 J _, V 
h m — - 

em 



ml 



dz 



l2^e 2 c a 2 C F a 2 s {l - r) 



mt 



7s + lnr - 21n(l 



MS 

2 



3m 2 s 2 



8(1 + r- z), 



x 



(16) 



where we have explicitly substituted the LO color-octet coefficient in (3) and the NLO 
correction to the color-octet matrix element (15). The subscript MS reminds that the color- 
singlet coefficient dF\ (//) is determined in accordance with the MS factorization scheme. 

Substituting the analytic expressions for the differential cc(}P^) production cross sec- 
tion, as assembled in (13), into (16), one can readily solve the desired differential color-singlet 
coefficient: 



dz 



64ne 2 a 2 C F aim, 



MS 

(1 -r)6(l + r-z) 



9iV> 2 




2z 



4r 



z — 2r 



1 + r 



+ 



z 2 — Ar 
z — 2r 



ml dcrfir 
6N C dz 



+ Vz 2 - 4r 



x 



(17) 



4 Note the authors of Ref. [9] used DR to regularize the UV divergence but a gluon mass to regulate the 
IR divergence for the perturbative color-octet NRQCD matrix element. Here we use DR to regulate both 
UV and IR divergences. 
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where daa n /dz is given in (13b). Although this short- distance coefficient is now free of IR 
singularity, as it should, it now depends logarithmically on the NRQCD factorization scale 
fi. However, when taking into account the /i-dependence of the color-octet matrix element 
(Og c ( 1 5'o)), the physical h c production rate in (1) in principle does not depend on this 
artificial scale. 



C. The integrated short-distance coefficients F X LH and Fg U 

Inspecting the differential coefficients in (4) and (17), one finds that the distribution 
of h c becomes singular near its maximum energy. However, it was noticed long ago that 
NRQCD expansion breaks down near the kinematic boundary of quarkonium momentum 
distribution [26], thus the NRQCD prediction at fixed order is no longer trustworthy. In 
order to reliably describe the energy distribution near the kinematic end point, one should 
incorporate the effect of nonperturbative shape function [27], as well as resum large Sudakov 
logarithm [28] , consequently the quarkonium spectrum near the endpoint will turn over and 
get smeared. However, including such refinement is beyond the scope of this work, which 
well deserves a dedicated study. 

On the other hand, the integrated h c production rate is much less sensitive to the above- 
mentioned effects, and NRQCD velocity expansion is believed to work well for this quantity. 
Moreover, the total production rate of the h c is also experimentally accessible. Therefore, it 
is of phenomenological incentive to find the integrated short-distance coefficients -F\ LH and 

pLH 
^8 • 

Integrating (17) over z may seem straightforward, but it turns out to be difficult to obtain 
the analytic expression for the part involving da& n /dz. We utilize a simple trick [29], i.e., 
restarting from (12b), but this time integrating over z first, then followed by integrating over 
x\. After some algebras, we end up with the following integrated short- distance coefficients: 



7 + 7r- 



6Ane 2 a 2 C F a 2 m 



9N r s 2 



'.1-r) 



-ln-£s + 21n(l-r)- 



Am 2 



65 - 84r 
12(1 -r) 



9r 2 



-I In r + 

6(1 -r) 2 + 

327r 2 e 2 a 2 a q m r 



r(5 - 7r) In 



2 1+VT- 



l- 



n- 



16(1 -r) 1 



+ 



(14-15r)ln± 



8(1 -r) 3 / 2 



3s 2 



:i-r) 



(18a) 
(18b) 



The color-octet coefficient F 8 LH can be trivially deduced by integrating (4) over z. 

It is natural to take the factorization scale /i around m c . We then find that F 1 LH (//) 
becomes negative in most of the allowed range of r (including r 0.08 of phenomenolog- 
ical interest), except in a narrow window where r is very small. This has an immediate 
consequence, that the color-octet channel becomes indispensable if one wishes to predict a 
positive total production rate for e + e~ — > h c + light hadrons. 

It is enlightening to examine the asymptotic behaviors of (18) in the limit ^> m: 



MS 



,-iLH 



asym 



asym 



&A^e 2 a 2 C F a 2 m c 

9N c s 2 
32n 2 e 2 a 2 a s m c 

3^2 ' 



— lnr 
12 



In 



Am 2 



65 7, n 
— + -ln2 
12 2 



(19a) 
(19b) 
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The total production rate for e + e _ — > h c + light hadrons scales with the center-of-mass 
energy as 1/s 2 , which is the same as that for e + e~ — > J/ip+ light hadrons [29]. Never- 
theless, it is interesting to note that the leading scaling violation in (19a) is represented 
by a single-logarithmic term (oc lnr), while that in the process e + e _ — > J/tp + g g is by a 
double-logarithmic term (oc ln 2 r) [29]. Since the newly proposed perturbative QCD factor- 
ization program for quarkonium production [30] is based on m 2 c /s expansion, it is natural to 
envisage that refactorizing the higher-twist two-parton fragmentation function in Ref. [30] 
may provide a natural framework to identify these logarithms at 0(a 2 ) and resum them to 
all orders in a s . 



IV. h c PRODUCTION ASSOCIATED WITH CHARMED HADRONS 

In this section, we investigate the inclusive h c production associated with the charmed 
hadrons at B factories, i.e., e + e~ — > h c + X c5 . Our central task is again to infer two short- 
distance coefficients appearing in the NRQCD factorization formula (1). To avoid confusion, 
we will always associate a superscript "Charm" to F\ and F 8 in this section, to distinguish 
from the analogous coefficients associated with the process e + e" — > h c + light hadrons in 
Sec. III. 

Both color-singlet and octet channels first occur at 0(a 2 ), characterized with the parton 
processes e + e~ — > cc^P^, 1 5'q 8 ' ) ) + cc. Since both channels share the common parton 
kinematics, we list some useful 3-body phase-space formulas here. 

We denote the momenta of the cc pair, c, and c by P, ki, k 2 , respectively, with P 2 4m 2 , 
k\ — k\ — m 2 . Analogous to (5), we also introduce three fractional energy variables z, x\ 
and x 2 , which are subject to the constraint X\ + x 2 + z = 2. 

Since no massless particles are involved in the final state, the parton cross sections in 
both channels do not exhibit any IR singularity. Unlike in Sec. Ill, we thus can perform 
the calculation directly in 4 dimensions. Consequently, suffice it to know the 4-dimensional 
3-body phase space measure for the process 7* — > cc(P) + c{ki) + c(k 2 ): 

d*s = ^ I dz [ 1 dx ± . (20) 



Jx 



The integration limits of z have been explicitly specified, while those for the fractional energy 
of c, xf, read 



A. Deducing dF^ ha 



We start with calculating the differential color-singlet coefficient dF^ harm . At the low- 
est order in a s , only four Feynman diagrams need to be considered for the process 7* — > 
cc^P^^ + cc. This calculation is quite similar to the analogous one for e + e~ — > J/tp + cc [31] 
and e + e~ — > Xcjijlc) + cc ( J = 0, 1,2) [32]. The perturbative matching calculation for this 
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coefficient is rather straightforward, so we directly present the result: 
6 l^ ;n-n: | 7(1 + r - z){\ - z){z 2 - 4r) 



dFf ha 
dz 



243m c s(2 - z)% 4 



768r 4 - 384r 3 (8 + 5r)z 



(2 - *) 4 

-64r(8 - 16r - 128r 2 - 35r 3 )^ 2 + 32r(56 - 64r - 310r 2 - 43r 3 )z 3 

+ 16(8 - 140r + 136r 2 + 438r 3 + 25r 4 )^ 4 - 8(48 - 80r + 136r 2 + 336r 3 - 15r 4 )/ 



+4(152 + 164r - 112r 2 + 132r 3 - 7r 4 ),2 e - (672 + 560r - 400r 2 + 84r 3 - 6r 4 )/ 



ii 



+2(300 + 182r - 60r 2 - r s )z* - 8(49 + 14r - 2r 2 )^ 9 + (130 + I7r)z w - 18z 

r , zJl + r - z + J(l - z)(z 2 - 4r) 
-— In — - y x 

2z Zy /i + r - z - ^(l - z){z 2 - 4r) 

- 192r 4 + 96r 3 (8 + 3r)z + 32r(4 - 8r - 39r 2 - 6r 3 )^ 2 - 32r(10 - 4r - 26r 2 - r 3 )z 3 
4(8 - 132r - 44r 2 + 48r 3 + 5r 4 )z 4 + 2(48 - 192r - 64r 2 + 48r 3 + 3r 4 )/ 



2(72 + 36r + 82r 2 + llr i )z b + 4(38 + 61r + 13r 2 y - (90 + 59r)z s + lQz 



(22) 



B. Deducing dF 8 Charm 

Next we proceed to calculate the differential color-octet coefficient <iF 8 Charm . At the order 
a 2 , we need consider six Feynman diagrams for the parton process 7* — > cc^Sy 1 ) + cc. The 
perturbative matching for this coefficient is analogous to that in Sec. Ill A, which is also 
quite straightforward. Here we just give the result: 



dR 



Charm 



1 



X 



2*^(1 -z){z 2 -4r) 

(1 : 27m c s (1 + r - z)(2 - z) 2 z 3 \ 3(2 - *) 4 (1 + r - zfl 2 

96r 3 (l + rf - 96r 2 (l + r) 2 (4 + 8r + r 2 )z 

+ 16(6 + 59r + 186r 2 + 328r 3 + 284r 4 + 51r 5 - 2r 6 )z 2 
-8(48 + 374r + 770r 2 + 965r 3 + 551r 4 + 41r 5 - r 6 )z 3 
+ (592 + 3392r + 5556r 2 + 5806r 3 + 2234r 4 + 90r 5 - 6r 6 )* 4 
-2(232 + 704r + 1162r 2 + 1131r 3 + 244r 4 - 5r 5 )* 5 

+ (262 - 207r + 217r 2 + 293r 3 + llr 4 )/ - (184 - 397r - 186r 2 - 15r 3 )z 7 



+ (116 - 139r - 40r 2 )z 8 - (40 - 13r)* 9 + 6z 
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, z\J\ + r - z + 7(1 - z)(z 2 - 4r) 
+ r In v . 7V ' x 



8r 3 (l + r) - 8r 2 (4 + 5r)z - 2(4 - 2r - 22r 2 - r 



z 7l + r - 2; - 7(1 - 2 )(^ 2 - 4r) 



+ (8 + 8r + 40r 2 - 6r 3 )* 3 - (14 + 33r - 5r 2 )z 4 + (2 - 19r)z 5 + I2z 



(23) 
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C. Fragmentation function for c — > h c and integrated cross section 



At first sight, the differential coefficients in (22) and (23) may look too disordered to 
extract anything useful. However, this is just a disguise, since we are certain that in the 
asymptotic limit y/s ^> m c , the differential h c production rate associated with cc must be 
dominated by the fragmentation mechanism: 

da[e+e - _> hc{P) + x cc ] 

— = 2aD Cr + hc {z), (24) 

where a = N c 7r ^ Q is the cross section for the process e + e~ — > cc, and D c ^ hr (z) stands for 
the fragmentation probability for c into the h c carrying the energy fraction z. The factor 2 
arises because both c and c can fragment into h c with equal probability. 

According to the NRQCD factorization ansatz, the fragmentation function of c into h c 
can be refactorized as 

D „ lc ( 2) = <.( 2 )M^ + ^(,)Ml^)), (25) 

!ll c !ll c 

where d!{ c (z) and d$ c (z) are the corresponding short- distance coefficient (jet) functions. 

Comparing (24) with (1), one can immediately realize that, these short-distance functions 
can be obtained by taking the asymptotic (r — > 0) limit of dF^ harm /dz (n = 1, 8): 



m, 



dFf ha 



~ CT dz nu 

16a 2 s z{l - z) 2 {64 - 128^ + 176z 2 - 160z 3 + UOz 4 - 56z 5 + 9z 6 



d h 8 <(z) 



dFp ha 



2& dz 



243(2 - zf 
a 2 s z(l - z) 2 (48 + 8z 2 - 8^ 3 + 3z 4 
asym " 162(2 - zf 



(26a) 
(26b) 



Both the expressions for d^ c (z) and d^ c (z) in (26) fully agree with Ref. [33]. 

Honestly speaking, the B factory energy is far from being asymptotically large. Therefore 
for this phenomenologically relevant case, the fragmentation function calculated in (26) can 
hardly faithfully reproduce the energy distribution of the h c depicted in (22) and (23). 

Like what has been done in Sec. Ill C, it is also of both theoretical and phenomenolog- 
ical interest in knowing the integrated production rate for e + e~ — > h c + charmed hadrons. 
Conceivably, it seems extremely challenging, if not impossible, to complete the integration 
of (22) and (23) over z in closed form. Nevertheless, it is quite easy to infer the asymptotic 
behavior of the total cross section with the help of the fragmentation function: 



a[e + e — >■ h c + X c . 



2d" / dzD c ^ hc (z) 
Jo 



16a 2 cr 



asym 

18107 - 26110 In 2 (^(^i)} 773 - 1110 In 2 (C>£ c (%)) 
8505 ^ c + 12960 mf 



(27) 



In contrast to (19), the production rate for h c + charmed hadron exhibits much slower asymp- 
totic decrease (oc 1/s) than that for h c + light hadrons (oc 1/s 2 ), which clearly corroborates 
the dominance of the fragmentation mechanism at high energy (pr)- 
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V. PHENOMENOLOGY 



With various color-singlet and color-octet short-distance coefficients determined in 
Sees. Ill and IV, we are ready to make a concrete analysis for the inclusive h c production 
at B factories and assess its observation prospects. 



A. Input parameters 

At the B factory energy, we choose the running QED coupling constant Oi(-\fs) = 1/131, 
and the running QCD coupling constant a s (y/s/2) 0.211. For the process e + e~ — > h c + 
light hadrons, we take the occurring factorization scale \i = m c . An important source of the 
uncertainty in our predictions is rooted in the uncertainty about charm quark mass. We 
take the customary value m c =1.5 GeV, but allow it to float between 1.3 to 1.8 GeV in 
order to closely assess this uncertainty. 

We need also specify the values of two nonperturbative NRQCD matrix elements appear- 
ing in the factorization formula (1). Upon vacuum saturation approximation, one can relate 
the color-singlet matrix element with the square of the first derivative of the radial wave 
function at the origin for IP charmonium, which is calculable in quark potential models [14]: 

(OtCPi)) « 2TPi {0ieJ(SPj)) * ^'^ (0)|2 - (28) 

If the value of R' 1P (0) is calculated from the Buchmuller-Tye potential model [34], one then 
finds ((^OPi)) = 0.322 GeV 5 . 

In contrast to the color-singlet matrix element, no reliable phenomenological model cal- 
culations are available for the color-octet matrix element (Og c ( 1 S'o)), since neither vacuum 
saturation approximation nor the quark potential models are applicable in this situation. 

Fortunately, the approximate heavy quark spin symmetry in NRQCD can be invoked to 
connect the color-octet matrix elements of h c and Xcj (J = 0, 1, 2) [14]: 

(O h 8 tS )) « ^—(0^(350). (29) 

There have been available a number of phenomenological studies for inclusive Xcj production 
in hadron collision or in B decay experiments [9, 35-39], and the color-octet matrix elements 
(Og CJ ( 3 S\)) have been fitted by various groups over years. We can use (29) to translate 
their fitted color-octet matrix elements for Xcj to the desired one for h c . In Table I, we have 
enumerated some values of (O^^So)) excerpted from various references. 

Alternatively, one can infer an order-of-magnitude estimate for the lower bound on 
(Og c ( 1 S'o)) using the renormalization group equation (RGE), which governs the renormal- 
ization scale dependence of this color-octet matrix element. The solution to the RGE at 
leading order in a s reads [14] 

(Os ( S )) mc = (O s ( So)), + ^ In j m , , (30) 

where /3 = UNc 3 2nf = 9 is the one-loop coefficient of QCD (3 function with rif — 3 light 
quark flavors. The subscript of the color-octet matrix element specifies the renormalization 
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TABLE I: Variation of the predicted total cross sections for e + e~ — > h c + light hadrons and for 
e + e _ —7- h c + charmed hadrons with the value of (Og c ( 1 5'o)) (in units of GeV 3 ). We have taken 
a = 1/131, a s (\/s/2) = 0.211, m c = 1.5 GeV, and the color-singlet matrix element = 
0.322 GeV 5 . The color-octet matrix element is taken from various references, while the last entry 
provides a lower bound inferred from the renormalization-group equation running in (31). 



Cross Section~~~~~~----~____^^ 


Ref. [35, 36] 
0.009 - 0.01 


Ref. [37] 
0.022 - 0.025 


Ref. [9, 38] 
0.014 - 0.020 


Ref. [39] 
0.039 


RGE 
> 0.0085 


e+e" -> he + X LU (fb) 


86.7- 97.6 


229.2 - 262.0 


141.5 - 207.2 


415.5 


> 81.2 


e+e" -> h c + X cE (fb) 


9.9-10.0 


10.2 - 10.3 


10.1 - 10.2 


10.7 


> 9.9 



e + e -* & c +light hadrons e + e -> /i e +charmed hadrons 




FIG. 1: The energy distribution of h c in the processes e + e — > h c + Vlh (left panel) and e + e — > 
h c + X c5 (right panel) at y/s = 10.58 GeV. The color-octet matrix element (Og c ( 1 S'o)) is taken as 
0.02 GeV 3 . Three curves in each plot correspond to taking m c = 1.3, 1.5, and 1.8 GeV, respectively. 



scale affiliated with the operator C 8 c ( 1 S'o). Taking jj, = m c v, and assuming the matrix 
element (0 8 c ( 1 So))m c v is nonnegative, one then gets: 



{0aCs „ )W> _^fM)wm. (31) 

243 V aJm r ) J mi 



c 



Taking a s (m c ) ~ 0.35 and a s (m c v) ~ v ~ 0.55, we then obtain ((9 8 ( 1 S'o)) mc > 0.0085 GeV 3 . 
This estimate is of course a very rough one. Nevertheless, as one can clearly see from Table I, 
all the phenomenologically determined values for the matrix element (Og c ( 1 S'o)) mc seem to 
be compatible with this bound. 

In the following numerical analyses, we will take (Og c ( 1 S'o)) mc = 0.02 GeV 3 , a medium 
value among those tabulated in Table I. 



B. Numerical results 

Substituting the differential short-distance coefficients dF^ R (n = 1,8) given in (4) and 
(17) into (1), we obtain the energy distribution of h c in the process e + e _ — > h c + light 
hadrons at B factory energy, as shown in the left panel of Fig. 1. Similarly, substituting the 
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differential coefficients dF^ harm (n = 1,8) given in (22) and (23) into (1), we then obtain 
the energy spectrum of h c in the process e + e~ — > h c + charmed hadrons, as shown in the 
right panel of Fig. 1. 

As can be clearly seen from Fig. 1, the h c energy spectra in two production channels are 
markedly different. In the former case, the differential production rate of h c sharply rises 
and diverges at the maximal energy of h c ; while in the latter, the distribution turns over 
and vanishes as the energy of h c approaches its maximum. As was discussed in Sec. Ill C, 
for the h c production associated with light hadrons, our prediction to the high-z part of 
the h c spectrum becomes untrustworthy, due to the breakdown of perturbative and velocity 
expansions near the endpoint region. An appropriate treatment requires resumming Sudakov 
logarithms as well as incorporating nonperturbative shape function, consequently one then 
expects the h c spectrum will be smeared, and turn over near the upper endpoint, rather 
than diverge. 

As we will see shortly, this process is largely dominated by the color-octet channel rather 
than the -singlet channel. From (4), the h c spectrum in color-octet channel at LO in a s is a 
sharp 5-function spiked on the maximal h c energy. Even though taking into account that the 
radiation of soft gluons would smear the h c energy spectrum in this channel, it is reasonable 
to expect that the majority of the h c events will still be located near the upper end of the 
momentum spectrum. This may serve as some useful guidance for the experimentalists. 

Our predictions for the integrated cross sections should be much more reliable than 
the differential distributions. Substituting the integrated short- distance coefficients 
(n = 1,8), which are collected in (18), into (1), we obtain the total h c production rate 
for the process e + e~ — > h c + light hadrons; for the channel e + e _ — > h c + charmed hadrons, 
the total cross section can be reached by numerically integrating the respective differential 
distribution. 

With our default value for the color-octet matrix element {Og c (*- So)) = 0.02 GeV 3 , we 
estimate the total cross section for h c + light hadrons to be about 207 fb, and that for h c + 
charmed hadrons to be about 10 fb. In Table I, we tabulate various predictions for the h c 
total cross sections in both production channels by adopting different values of color-octet 
matrix elements, which are in the following ranges: 

a[e + e- ^h c + X LH ] = 81.2 - 415.5 fb, (32a) 
a[e + e~ ^h c + X cc -] = 9.9 - 10.7 fb. (32b) 

One immediately observe that, the cross section of the former channel is quite sensitive 
to the value of (Og c ( 1 S'o)), but that of the latter is not sensitive to it at all. This clearly 
indicates that the former process is dominated by the color-octet channel, while the latter 
is dominated by the color-singlet channel. 

It is interesting to contrast the inclusive h c production at B factories with the anal- 
ogous inclusive J/tp production processes, which have been recently measured by Belle 
Collaboration [15]: 

a[e + e- ->■ J/ip + X LH ] = 430 ± 90 ± 90 fb, (33a) 
a[e + e~ -> J/rj) + X c5 ] = 740 ± 80±|g fb. (33b) 

We see that the cross section for h c + light hadrons is comparable in magnitude with the 
observed production rate for J/ip+ light hadrons. In this work we have only implemented 
the color-octet short- distance coefficient F 8 LH at LO in a s . It was recently discovered that 
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e + e -> /i c +light hadrons e + e -» /i c +charmed hadrons 




m e (GeV) m c (GeV) 



FIG. 2: The total production rate for h c as a function of m c at y/s = 10.58 GeV. The left panel is 
for the process e + e _ — > h c + Vlh.) and the right panel is for e + e~ — > h c + X c g. The band in each 
plot characterizes the uncertainty estimated by varying the color-octet matrix element (O s c ( l So)) 
from 0.0085 to 0.039 GeV 3 , where the central curve corresponds to fixing (O^^Sq)) at 0.02 GeV 3 . 



there may exist a large positive 0(a s ) correction to this coefficient [40]. If we include this 
perturbative correction, the total production rate for h c + light hadrons may easily reach 
the value given in (33a). Therefore, copious h c + light hadrons events should already have 
been produced at B factories. 

One naturally expects that the NRQCD factorization framework will eventually break 
down as one keeps pushing down the center-of-mass energy. Leaving this caveat aside, it may 
still be tempting to bluntly apply our formula to the process e + e~ — > /i c 7r + 7r~ at y/s = 4.17 
GeV, first observed in the CLEO experiment [7]. Taking the default value of the color-octet 
matrix element (£>g c (%)) = 0.02 GeV 3 , and choosing a s (y/s/2) = 0.305, m c = 1.5 GeV, we 
then obtain <j[e + e~ — > h c + XlhJ ~ 2.5 pb, which is considerably lower than the measured 
cross section 15.6 ± 2.3 ± 1.9 ±3.0 pb [7]. This may indicate that, when y/s gets close to the 
open charm threshold, y/s ~ 2m c , the h c production is almost saturated by the exclusive 
events, therefore the NRQCD factorization, which is tailor-made to tackle the inclusive 
quarkonium production, becomes inevitably untrustworthy. To handle this situation more 
appropriately, one likely needs appeal to the even lower-energy effective field theory such as 
the potential NRQCD. 

In sharp contrast with inclusive J/ip production, our predicted production rate for h c + 
charmed hadrons at the B factories is about one order of magnitude smaller than that for 
h c + light hadrons. This is a normal hierarchy pattern, being consistent with the heuristic 
expectation based on the kinematics consideration. Practically speaking, the low production 
rate may render the experimental measurements of this production channel of h c difficult. 

One should caution that, the predicted h c production rates given in (32) are still subject 
to large theoretical uncertainties. Besides the value of (Og c (* S )} , the charm quark mass 
constitutes one major source of the uncertainties. In Fig. 2 we explicitly show how the 
total production rates for h c vary with m c . A useful message conveyed by Fig. 2 is that, 
the production rate for e + e~ — > h c + charmed hadrons is much more sensitive to m c than 
that for e + e~ — > h c + light hadrons (this point can also be seen in Fig. 1). This fact can 
be best explained by going to the asymptotic limit y/s ^> m c . In such a limit, the process 
e + e~ — > h c + X c5 is dominated by the color-singlet fragmentation mechanism, thus in light 
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e + e -> /i c +light hadrons e + e -» /i c +charmed hadrons 




V«" (GeV) (GeV) 



FIG. 3: The integrated /i c production cross section as a function of \/s for the processes e + e _ — >• 
/i c + Xlh (left panel) and e + e~ — >• /i c + X c£ (right panel). The band in each plot measures the 
uncertainty brought in by varying the color-octet matrix element (O g c ( l So)) from 0.0085 to 0.039 
GeV 3 , where the central curve corresponds to fixing (O^^So)) at 0.02 GeV 3 . 

of (27), one expects a oc (Oj lc ( 1 Pi))/m^. The process e + e~ — > h c + X^u is dominated by the 
color-octet channel. From (19b), one finds that, asymptotically a oc (Os( 1 5'i))/m c . It is this 
very different power-law scaling behavior that accounts for the much stronger sensitivity of 
a[e + e~ — > h c + A cE ] to charm quark mass. 

Finally, in Fig. 3 we show how the integrated h c production rates vary with the center- 
of-mass energy. As yfs increases, the total cross section of e + e~ — > h c + Alh descends much 
steeper than that of e + e _ — > h c + A cg . This can be easily attributed to the fact that at 
large s, the production rate in the former process scales as 1/s 2 , while that in the latter case 
scales only as 1/s, as dictated by the fragmentation mechanism. 

C. Observation prospects of inclusive h c production at B factories 

Thus far, there is yet no any experiment at B factories dedicated to measure the inclusive 
h c production rate. Perhaps it is partly due to the lack of phenomenological incentive, and 
more importantly, due to the difficulty of reconstructing the h c signals. 

Our prediction in (32), indicates that a large number of h c + light hadrons events should 
have already been produced at B factories, given the large integrated luminosity that has 
been accumulated by two B factory experiments near the T(4S) resonance. As we have 
shown before, the prospective measurements of this channel would provide a promising 
window to test the color-octet mechanism in quarkonium production, and impose some 
useful constraint on the color-octet matrix element (O^^Sq)}. In view of this, there should 
be sufficient phenomenological impetus for experimentalists to pursue the measurements. 

We can examine the observation potential for inclusive h c production more quantitatively. 
Up to present, Belle experiment has accumulated about 1000 ftT 1 data near the T(45*) 
resonance. Taking cr[e + e~ — > h c + Alh] ~ 100 — 400 fb from (32a), we thus estimate that 
roughly (1 — 4) x 10 5 h c events have been produced at Belle at y/s = 10.58 GeV. 

The two known decay channels of the h c meson are h c — > 2(7r + 7r~)7r° and h c — > 7^7. The 
multi-pion decay seems to be a clean and potentially useful tagging mode fr the h c signal. 
Nevertheless, due to the greater branching fraction of the latter channel, it might be of some 
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advantages by utilizing the El transition h c — » r? c 7 to reconstruct the h c signal, i.e., first 
reconstruct an i] c , then check if there exists a narrow peak around 3525 MeV in the 7 + i] c 
invariant mass distribution. 

The main experimental difficulty hinges on how to efficiently reconstruct the rj c . It seems 
to be a standard practice for the Belle Collaboration to reconstruct the r] c meson via several 
hadronic decay modes, e.g. K s K + n~ +c.c, n + it~ K + K~ , 2(K + K~), 2(7r + 7r~), 3(7r + 7r~) by 
looking for the charged tracks (In fact, they have employed this technique in searching for 
the radiative decay processes T(1S,2S) — >■ 7^ c [41, 42]). A conservative estimate gives that 
about the 1% r\ c events can be reconstructed 5 . 

Taking B(h c -> 777J « 50% [5], we then expect roughly (1 - 4) x 10 5 x 50% x 1% = 
500 — 2000 reconstructed h c events. 

A large number of reconstructed h c events seems to indicate a bright observation prospect, 
however, one must be alert to the potentially huge combinatorial background, since there 
are a lot of pions, kaons and photons in the events. A careful study of the background level 
is crucial from the experimental perspective. 

VI. SUMMARY 

In this work, we have studied the inclusive production of the h c meson associated with 
the light hadrons and with the charmed hadrons at the B factories, respectively, in the 
framework of NRQCD factorization. We have only considered the lowest-order contribution 
in v, which involves both the color-singlet channel and the color-octet channel. 

We have explicitly verified that, for e + e~ — > h c + light hadrons, including the color-octet 
channel is pivotal in rendering the IR-finite prediction for the inclusive h c production rate. 
Moreover, at B factory energy yfs = 10.58 GeV, within some reasonable choices for the 
color-octet matrix element (Og c ( 1 S'o)), we find that the total h c production rate varies in 
the range between 0.08 and 0.4 pb. Furthermore, this process is found to be dominated by 
the color-octet channel. Future measurement of this process at B factories may provide an 
explicit test on the color-octet mechanism, as well as put a useful constraint on the value of 
the color-octet matrix element (O^^Sq)). 

Since this process is dominated by the color-octet channel, it is reasonable to expect that 
the bulk of the h c events will populate near the upper end of the momentum spectrum. For 
the experimentalists to unambiguously measure the h c production rate at the B factories, it 
is perhaps imperative to first have a reliable prediction for the momentum spectrum of the 
h c 6 . This direction is certainly worth further investigations. 

For the process e + e _ — > h c + charmed hadrons, we find that the total production rate is 
almost saturated by the color-singlet channel. The predicted total production cross section 
at y/s = 10.58 GeV is about 10 fb if the charm quark mass is taken to be 1.5 GeV. This is 
about one order-of- magnitude smaller than the cross section for e + e~ — > h c + light hadrons. 



5 This estimate is obtained from the product of detection efficiency and the branching fractions of afore- 
mentioned rjc decay channels. With the efficiency taken with around 25%, it is reasonable to assume that 
1% r\ c events will be reconstructed. We thank C. Z. Yuan for providing this estimate and for explaining 
to us some experimental details. 

6 We thank C.-P. Shcn for stressing this point to us, and for having performed a detailed Monte Carlo study 
on reconstructing the h c events at the Belle experiment. 
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We note that this prediction bears large uncertainty, for it is quite sensitive to the charm 
quark mass. If allowing m c to float from 1.8 GeV to 1.3 GeV, this cross section would vary 
from a few fb to 30 fb. It is interesting to see whether the future experiments can observe 
this process or not. 

Note added. After this paper was submitted, there has recently appeared a related work in 
arXiv [46], which also studied the inclusive h c production in e + e~ annihilation in NRQCD 
factorization framework. For the h c production associated with the light hadrons, these 
authors adopted a factorization scheme different from the MS scheme, and calculated the 
integrated NRQCD short-distance coefficients using some numerical recipe. 
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Appendix A: Isolating the infrared divergence for e + e~ — > ccQ-P^) + gg in dimen- 
sional regularization 

In this Appendix, we explain how to calculate the IR-divergent integral in (12a) in di- 
mensional regularization. The D- dimensional 3-body phase space measure d& 3 and the 
integrand Idiv(xi,z) have already been presented in (6) and (10), respectively. Inspecting 
these expressions, one can readily recognize that the IR singularity would arise from two 
distinct phase space corners: z — > 1 + r, x\ — > and z — > 1 + r, X\ — > 1 — r {x 2 — > 0), by 
integrating the first and second terms in (10), respectively. 

First note that both the integrand and the phase space measure in (12a) are symmetric 
under the interchange x\ -H- x 2 . A useful shortcut to perform the Xi-integral is to integrate 
over only the lower half of its allowed range, then multiply the result by 2: 

na(z)+b(z) pa(z) 

/ dx 1 = 2 dx 1} (Al) 

J a(z)—b(z) J a(z)—b(z) 

where a(z) and b(z) are defined in (8). As a simplification, the IR singularity now is solely 
caused by the first term in (10), when integrated over the phase space corner z — > 1 + r, x 1 — > 
0. The integral involving the second term now becomes IR finite, thereby can be worked 
ont directly in 4 dimensions. 

We then face the IR-divergent integral of the following type: 

Hz) - r (Z) 9{XI > Z) (A2) 

{) ~ Ja^-Kz) dX \l + r-z- Xl fxf{l - cos2 0Y ' (A2j 

where cos9 as a function of z and x\ has been given in (7), and g(xi,z) is an arbitrary 
function that is regular at z — 1+r. Note that the e-dependent factors in the integrand come 
from the D-dimensional 3-body phase space measure, playing the role of the IR regulator. 

A brute-force calculation of I(z) while keeping the full e dependence is a challenging task, 
if not impossible. It is desirable if one can make the IR divergence explicit prior to carrying 
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out the x\ integration. This is indeed feasible, and in the analysis of the decay process 
XbJ ccg, Bodwin et al. [25] have elaborated on how to fulfill such a goal by utilizing a 
simple trick. In below, we will employ a strategy that is closely analogous to theirs 1 . 

The key is to realize that x\x\ = a 2 (z) — b 2 (z) = l+r — z. As suggested by this identity, 
it may seem advantageous to introduce a new integration variable t: 



t = 



1 + r 



Xi 



= a(z) + b(z) cos 9. 



(A3) 



The integration range for this new variable turns out to be t < t < a(z) + b(z), where 

2(1 + r- z) 



(A4) 



In term of the new variable t, equation (A2) becomes 



I(z) = 



1 



[1 + r - z) 



l+2e 



a(z)+b(z) 



,2e 



dt- 



to 



l-ty (i-cos 2 ey 



(A5) 



In light of (A3), cos# now should be understood as a function of t, a(z) and b(z). 

With the aid of the familiar identity about distributions, we can express (1 + r — z)~ 1 ~ 2t 

as 



1 



5(1 



[1 + r - z) 



l+2e 



2e (1 - y/F) 



Ac 



+ 



1 


-2e 


1 + r — z 


+ 



ln(l 



1 + r 



(A6) 



As advertised, we have successfully separated the IR pole that is accompanied with the 5(1 + 
r—z) function. The IR-regular remainders are partially encoded in the "+" -functions, whose 
integration property has been specified in (14). To the desired accuracy, the logarithmic "+" 
function is not needed in this work. 

One can be readily convinced that the integral in (A5) is finite in the limit e — > 0. 
Therefore, to our intended accuracy, it can be evaluated by simply Taylor-expanding the 
integrand through the first order in e. 

Using the master formulas (A5) and (A6), it is now a straightforward exercise to repro- 
duce the analytic expression of dtr div /dz as given in (13a). One can further integrate this 
expression over z to arrive at 



o"di 



128ire 2 a 2 C F a 2 (l - r) c £ (4tt) 



Sm 2 i 
+41n(l - r) - 



3r 



r(i-e) 

lnr — 1 



1 

em 



2 In 



Am 2 



(AT) 



Wc note that the process rib — > XcJ + 99 (J = 0, 1, 2) has recently been analyzed in NRQCD factoriza- 
tion [43], which shares essentially the identical kinematics as our process. The methodology of isolating 
the IR singularity in DR in [43] is similar to what is adopted in [25]. 
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Appendix B: r) c production associated with charmed hadrons in e+e annihilation 



In this Appendix, we consider e + e~ — > r) c + charmed hadrons at lowest order in v and a s 
(note that e + e~ — >■ r] c + gg is forbidden by the charge conjugation invariance). The result 
presented in this section can be viewed as a byproduct of the analysis made in Sec. IV. We 
note that, the process e + e~ — >■ i] c + cc has already been investigated some time ago [32], and 
our calculation may serve as an independent check. 

At the lowest order in v, the NRQCD factorization formula for the r\ c inclusive production 
reads [14]: 



da[e + e ->■ i] c + X cS \ = 



dFf ha 



mi 



-{OrCS )) + O(av 2 ), 



(Bl) 



where the color-singlet production operator Of^So) has been defined in [14], <iF 1 Charm is the 
corresponding short- distance coefficient. Unlike the P-wave charmonium production, here 
we are justified to ignore the color-octet channel since its relative importance is suppressed 
by v A . 

The perturbative matching calculation for the color-singlet short- distance coefficient is 
completely analogous to Sec. IV, and we directly present the result: 



dF^ ha 



32ire 2 c a 2 a 2 s 



2z 



[l-z)(z 2 -Ar) 



x 



dz 81m c s(2- z) 2 z 3 1 3(2 - zf V 1 + r-z 

96r 3 (l + r) - 96r 2 (4 + 6r + r 2 )z + 16(6 - 25r + 38r 2 + 43r 3 - 2r A )z 2 
-8(24 - 138r + 58r 2 + 51r 3 - r A )z 3 + (112 - 944r + 372r 2 + 118r 3 - 6r 4 )z 4 
-2(24 - 120r + 62r 2 + r 3 )z 5 + (54 + 5r + 13r 2 )^ 6 - (28 + llr)z 7 + 6^ 

+rln ( z ^ 1 + r ~ z + v / ( 1 ~ z)(z 2 -Ar)' 
\z^l + r - z - a/(1 - z)(z 2 - Ar) 

-4r(l +r)z 3 + (10 + r)^ 4 



8r d 



32r 2 z - 2(4 - 6r - r 3 )z 2 



(B2) 



Our expression seems to be considerably compact than its counterpart in [32]. We have 
numerically checked that once using their phenomenological input parameters, we can re- 
produce their predicted production rate for r\ c + X c5 at B factory energy. 

Similar to the fragmentation function of c into h c in (25), the fragmentation function of 
c into Tj c can also be factorized as 



rrvi 



(B3) 



where is the color-singlet coefficient function, which agrees with the one give in [45]. 

Following the steps described in Sec. IV C, it is straightforward to identify d T [ c (z) by 
taking the asymptotic limit to (B2): 



d?(z) = 



Charm 



m c dF[ 
2a dz 



asym 



16a 2 ^(l - zf (48 + 8^ 2 - 8^ 3 + 3z 4 ) 
243(2 - zf : 



(B4) 
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which fully agrees exactly with the fragmentation function of c — > i] c first obtained in [44, 45]. 
Obviously, this expression only differs with dg c in (26) by a color factor. 

Finally we assess the integrated cross section for e + e~ — > r] c + X c5 . It is difficult to obtain 
the analytic result by integrating (B2) over z. Nevertheless we are content with knowing its 
asymptotic behavior by resorting to fragmentation approximation: 



This expression is compatible with the one give in [45]. 
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